home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
QRZ! Ham Radio 8
/
QRZ Ham Radio Callsign Database - Volume 8.iso
/
pc
/
files
/
ant_nec
/
nec81tar.z
/
nec81tar
/
eksc.f
< prev
next >
Wrap
Text File
|
1991-05-13
|
14KB
|
579 lines
C $TITLE: 'EKSC'
C $NOFLOATCALLS
C
SUBROUTINE EKSC (S,Z,RH,XK,IJ,EZS,ERS,EZC,ERC,EZK,ERK)
C COMPUTE E FIELD OF SINE, COSINE, AND CONSTANT CURRENT FILAMENTS BY
C THIN WIRE APPROXIMATION.
REAL*8 CONX(2),XK,RKB2,RHK,SHK,SS,CS
COMPLEX*16 CON,GZ1,GZ2,GP1,GP2,GZP1,GZP2,EZS,ERS,EZC,ERC,EZK,ERK
COMMON /TMI/ ZPK,RKB2,IJX
EQUIVALENCE (CONX,CON)
DATA CONX/0.,4.771341189D0/
C**
IJX=IJ
ZPK=XK*Z
RHK=XK*RH
RKB2=RHK*RHK
SH=.5*S
SHK=XK*SH
SS=DSIN(SHK)
CS=DCOS(SHK)
Z2=SH-Z
Z1=-(SH+Z)
C**
CALL GX (Z1,RH,XK,GZ1,GP1)
CALL GX (Z2,RH,XK,GZ2,GP2)
C**
GZP1=GP1*Z1
GZP2=GP2*Z2
EZS=CON*((GZ2-GZ1)*CS*XK-(GZP2+GZP1)*SS)
EZC=-CON*((GZ2+GZ1)*SS*XK+(GZP2-GZP1)*CS)
ERK=CON*(GP2-GP1)*RH
CALL INTX (-SHK,SHK,RHK,IJ,CINT,SINT)
EZK=-CON*(GZP2-GZP1+XK*XK*CMPLX(CINT,-SINT))
GZP1=GZP1*Z1
GZP2=GZP2*Z2
IF (RH.LT.1.E-10) GO TO 1
ERS=-CON*((GZP2+GZP1+GZ2+GZ1)*SS-(Z2*GZ2-Z1*GZ1)*CS*XK)/RH
ERC=-CON*((GZP2-GZP1+GZ2-GZ1)*CS+(Z2*GZ2+Z1*GZ1)*SS*XK)/RH
C**
RETURN
1 ERS=(0.,0.)
ERC=(0.,0.)
RETURN
END
C
C
C
SUBROUTINE EKSCX(BX,S,Z,RHX,XK,IJ,INX1,INX2,EZS,ERS,EZC,ERC,EZK,
1ERK)
C COMPUTE E FIELD OF SINE, COSINE, AND CONSTANT CURRENT FILAMENTS BY
C EXTENDED THIN WIRE APPROXIMATION.
INTEGER*4 INX1,INX2
REAL*8 CONX(2),XK,RKB2,RHK,SHK,SS,CS
COMPLEX*16 CON,EZS,ERS,EZC,ERC,EZK,ERK,GZ1,GZ2,GRK1,GRK2
COMPLEX*16 GZP1,GZP2,GR1,GR2,GRP1,GRP2,GZZ1,GZZ2
COMMON /TMI/ ZPK,RKB2,IJX
EQUIVALENCE (CONX,CON)
DATA CONX/0.,4.771341189D0/
C**
IF (RHX.LT.BX) GO TO 1
RH=RHX
B=BX
IRA=0
GO TO 2
1 RH=BX
B=RHX
IRA=1
2 SH=.5*S
IJX=IJ
ZPK=XK*Z
RHK=XK*RH
RKB2=RHK*RHK
SHK=XK*SH
SS=DSIN(1.D0*SHK)
CS=DCOS(1.D0*SHK)
Z2=SH-Z
Z1=-(SH+Z)
A2=B*B
IF (INX1.EQ.2) GO TO 3
CALL GXX (Z1,RH,B,A2,XK,IRA,GZ1,GZP1,GR1,GRP1,GRK1,GZZ1)
GO TO 4
3 CALL GX (Z1,RHX,XK,GZ1,GRK1)
GZP1=GRK1*Z1
GR1=GZ1/RHX
GRP1=GZP1/RHX
GRK1=GRK1*RHX
GZZ1=(0.,0.)
4 IF (INX2.EQ.2) GO TO 5
CALL GXX (Z2,RH,B,A2,XK,IRA,GZ2,GZP2,GR2,GRP2,GRK2,GZZ2)
GO TO 6
5 CALL GX (Z2,RHX,XK,GZ2,GRK2)
GZP2=GRK2*Z2
GR2=GZ2/RHX
GRP2=GZP2/RHX
GRK2=GRK2*RHX
GZZ2=(0.,0.)
6 EZS=CON*((GZ2-GZ1)*CS*XK-(GZP2+GZP1)*SS)
EZC=-CON*((GZ2+GZ1)*SS*XK+(GZP2-GZP1)*CS)
ERS=-CON*((Z2*GRP2+Z1*GRP1+GR2+GR1)*SS-(Z2*GR2-Z1*GR1)*CS*XK)
ERC=-CON*((Z2*GRP2-Z1*GRP1+GR2-GR1)*CS+(Z2*GR2+Z1*GR1)*SS*XK)
ERK=CON*(GRK2-GRK1)
CALL INTX (-SHK,SHK,RHK,IJ,CINT,SINT)
BK=B*XK
BK2=BK*BK*.25
EZK=-CON*(GZP2-GZP1+XK*XK*(1.-BK2)*CMPLX(CINT,-SINT)-BK2*(GZZ2
1 -GZZ1))
RETURN
END
C
C
C
SUBROUTINE ROM2 (A,B,SUM,DMIN)
C
C FOR THE SOMMERFELD GROUND OPTION, ROM2 INTEGRATES OVER THE SOURCE
C SEGMENT TO OBTAIN THE TOTAL FIELD DUE TO GROUND. THE METHOD OF
C VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION IS USED. THERE ARE 9
C FIELD COMPONENTS - THE X, Y, AND Z COMPONENTS DUE TO CONSTANT,
C SINE, AND COSINE CURRENT DISTRIBUTIONS.
C
REAL*8 A,B,TMAG1,TMAG2,TR,TI,DMIN,Z,DZ,DZOT,S,EP,ZE,ZEND
COMPLEX*16 SUM,G1,G2,G3,G4,G5,T00,T01,T10,T02,T11,T20
INTEGER*4 NM,NS
DIMENSION SUM(9),G1(9),G2(9),G3(9),G4(9),G5(9),T01(9),T10(9),
1T20(9)
DATA NM,NTS,NX,N/65536,4,1,9/,RX/1.E-4/
C**
Z=A
ZE=B
S=B-A
IF (S.GE.0.) GO TO 1
WRITE(*,18)
STOP
1 EP=S/(1.D4*NM)
ZEND=ZE-EP
DO 2 I=1,N
2 SUM(I)=(0.,0.)
NS=NX
NT=0
CALL SFLDS (Z,G1)
3 DZ=S/NS
IF (Z+DZ.LE.ZE) GO TO 4
DZ=ZE-Z
IF (DZ.LE.EP) GO TO 17
4 DZOT=DZ*.5
CALL SFLDS (Z+DZOT,G3)
CALL SFLDS (Z+DZ,G5)
5 TMAG1=0.
TMAG2=0.
C
C EVALUATE 3 POINT ROMBERG RESULT AND TEST CONVERGENCE.
C
DO 6 I=1,N
T00=(G1(I)+G5(I))*DZOT
T01(I)=(T00+DZ*G3(I))*.5
T10(I)=(4.*T01(I)-T00)/3.
IF (I.GT.3) GO TO 6
TR=REAL(T01(I))
C TI=AIMAG(T01(I))
TI=IMAG(T01(I))
TMAG1=TMAG1+TR*TR+TI*TI
TR=REAL(T10(I))
C TI=AIMAG(T10(I))
TI=IMAG(T10(I))
TMAG2=TMAG2+TR*TR+TI*TI
6 CONTINUE
TMAG1=DSQRT(TMAG1)
TMAG2=DSQRT(TMAG2)
CALL TEST(TMAG1,TMAG2,TR,0.D0,0.D0,TI,DMIN)
IF(TR.GT.RX)GO TO 8
DO 7 I=1,N
7 SUM(I)=SUM(I)+T10(I)
NT=NT+2
GO TO 12
8 CALL SFLDS (Z+DZ*.25,G2)
CALL SFLDS (Z+DZ*.75,G4)
TMAG1=0.
TMAG2=0.
C
C EVALUATE 5 POINT ROMBERG RESULT AND TEST CONVERGENCE.
C
DO 9 I=1,N
T02=(T01(I)+DZOT*(G2(I)+G4(I)))*.5
T11=(4.*T02-T01(I))/3.
T20(I)=(16.*T11-T10(I))/15.
IF (I.GT.3) GO TO 9
TR=REAL(T11)
C TI=AIMAG(T11)
TI=IMAG(T11)
TMAG1=TMAG1+TR*TR+TI*TI
TR=REAL(T20(I))
C TI=AIMAG(T20(I))
TI=IMAG(T20(I))
TMAG2=TMAG2+TR*TR+TI*TI
9 CONTINUE
TMAG1=DSQRT(TMAG1)
TMAG2=DSQRT(TMAG2)
CALL TEST(TMAG1,TMAG2,TR,0.D0,0.D0,TI,DMIN)
IF(TR.GT.RX)GO TO 14
10 DO 11 I=1,N
11 SUM(I)=SUM(I)+T20(I)
NT=NT+1
12 Z=Z+DZ
IF (Z.GT.ZEND) GO TO 17
DO 13 I=1,N
13 G1(I)=G5(I)
IF (NT.LT.NTS.OR.NS.LE.NX) GO TO 3
NS=NS/2
NT=1
GO TO 3
14 NT=0
IF (NS.LT.NM) GO TO 15
WRITE(*,19) Z
GO TO 10
15 NS=NS*2
DZ=S/NS
DZOT=DZ*.5
DO 16 I=1,N
G5(I)=G3(I)
16 G3(I)=G2(I)
GO TO 5
17 CONTINUE
C**
RETURN
C
18 FORMAT (' ERROR - B LESS THAN A IN ROM2')
19 FORMAT (33H ROM2 -- STEP SIZE LIMITED AT Z =,1P,E12.5)
END
C
C
C
SUBROUTINE GX (ZZ,RH,XK,GZ,GZP)
C SEGMENT END CONTRIBUTIONS FOR THIN WIRE APPROX.
REAL*8 XK,R,R2,RK,CRK,SRK
COMPLEX*16 GZ,GZP
C**
R2=ZZ*ZZ+RH*RH
R=DSQRT(R2)
RK=XK*R
CRK= DCOS(RK)
SRK=-DSIN(RK)
GZ=DCMPLX(CRK,SRK)/R
GZP=-DCMPLX(1.,RK)*GZ/R2
END
C
C
C
SUBROUTINE GXX (ZZ,RH,A,A2,XK,IRA,G1,G1P,G2,G2P,G3,GZP)
C SEGMENT END CONTRIBUTIONS FOR EXT. THIN WIRE APPROX.
REAL*8 XK,R,R2,R4,RK,RK2,RH2
COMPLEX*16 GZ,C1,C2,C3,G1,G1P,G2,G2P,G3,GZP
R2=ZZ*ZZ+RH*RH
R=DSQRT(R2)
R4=R2*R2
RK=XK*R
RK2=RK*RK
RH2=RH*RH
T1=.25*A2*RH2/R4
T2=.5*A2/R2
C1=CMPLX(1.,RK)
C2=3.*C1-RK2
C3=CMPLX(6.,RK)*RK2-15.*C1
GZ=CMPLX(DCOS(RK),-DSIN(RK))/R
G2=GZ*(1.+T1*C2)
G1=G2-T2*C1*GZ
GZ=GZ/R2
G2P=GZ*(T1*C3-C1)
GZP=T2*C2*GZ
G3=G2P+GZP
G1P=G3*ZZ
IF (IRA.EQ.1) GO TO 2
G3=(G3+GZP)*RH
GZP=-ZZ*C1*GZ
IF (RH.GT.1.E-10) GO TO 1
G2=0.
G2P=0.
RETURN
1 G2=G2/RH
G2P=G2P*ZZ/RH
RETURN
2 T2=.5*A
G2=-T2*C1*GZ
G2P=T2*GZ*C2/R2
G3=RH2*G2P-A*GZ*C1
G2P=G2P*ZZ
GZP=-ZZ*C1*GZ
RETURN
END
C
C
C
SUBROUTINE INTX (EL1,EL2,B,IJ,SGR,SGI)
C
C INTX PERFORMS NUMERICAL INTEGRATION OF EXP(JKR)/R BY THE METHOD OF
C VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION. THE INTEGRAND VALUE
C IS SUPPLIED BY SUBROUTINE GF.
C
INTEGER*4 NM,NS
REAL*8 T01R,T10R,TE1R,T01I,T10I,TE1I,T11R,T20R,TE2R,TE2I,T11I,
1 T20I,TE21,T00R,T00I,G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,
2 RKB2,B,S,EL1,EL2,Z,ZE,ZP,DZ,DZOT
COMMON /TMI/ ZPK,RKB2,IJX
DATA NX,NM,NTS,RX/1,65536,4,1.E-4/
C**
Z=EL1
ZE=EL2
IF (IJ.EQ.0) ZE=0.
S=ZE-Z
FNM=NM
EP=S/(10.*FNM)
ZEND=ZE-EP
SGR=0.
SGI=0.
NS=NX
NT=0
C**
CALL GF (Z,G1R,G1I)
C**
1 FNS=NS
DZ=S/FNS
ZP=Z+DZ
IF (ZP-ZE) 3,3,2
2 DZ=ZE-Z
IF (ABS(DZ)-EP) 17,17,3
3 DZOT=DZ*.5
ZP=Z+DZOT
CALL GF (ZP,G3R,G3I)
ZP=Z+DZ
CALL GF (ZP,G5R,G5I)
4 T00R=(G1R+G5R)*DZOT
T00I=(G1I+G5I)*DZOT
T01R=(T00R+DZ*G3R)*0.5
T01I=(T00I+DZ*G3I)*0.5
T10R=(4.0*T01R-T00R)/3.0
T10I=(4.0*T01I-T00I)/3.0
C
C TEST CONVERGENCE OF 3 POINT ROMBERG RESULT.
C
CALL TEST (T01R,T10R,TE1R,T01I,T10I,TE1I,0.D0)
IF (TE1I-RX) 5,5,6
5 IF (TE1R-RX) 8,8,6
6 ZP=Z+DZ*0.25
CALL GF (ZP,G2R,G2I)
ZP=Z+DZ*0.75
CALL GF (ZP,G4R,G4I)
T02R=(T01R+DZOT*(G2R+G4R))*0.5
T02I=(T01I+DZOT*(G2I+G4I))*0.5
T11R=(4.0*T02R-T01R)/3.0
T11I=(4.0*T02I-T01I)/3.0
T20R=(16.0*T11R-T10R)/15.0
T20I=(16.0*T11I-T10I)/15.0
C
C TEST CONVERGENCE OF 5 POINT ROMBERG RESULT.
C
CALL TEST (T11R,T20R,TE2R,T11I,T20I,TE2I,0.D0)
IF (TE2I-RX) 7,7,14
7 IF (TE2R-RX) 9,9,14
8 SGR=SGR+T10R
SGI=SGI+T10I
NT=NT+2
GO TO 10
9 SGR=SGR+T20R
SGI=SGI+T20I
NT=NT+1
10 Z=Z+DZ
IF (Z-ZEND) 11,17,17
11 G1R=G5R
G1I=G5I
IF (NT-NTS) 1,12,12
12 IF (NS-NX) 1,1,13
C
C DOUBLE STEP SIZE
C
13 NS=NS/2
NT=1
GO TO 1
14 NT=0
IF (NS-NM) 16,15,15
15 WRITE(*,20) Z
GO TO 9
C
C HALVE STEP SIZE
C
16 NS=NS*2
FNS=NS
DZ=S/FNS
DZOT=DZ*0.5
G5R=G3R
G5I=G3I
G3R=G2R
G3I=G2I
GO TO 4
17 CONTINUE
IF (IJ) 19,18,19
C
C ADD CONTRIBUTION OF NEAR SINGULARITY FOR DIAGONAL TERM
C
18 SGR=2.*(SGR+DLOG((DSQRT(B*B+S*S)+S)/B))
SGI=2.*SGI
19 CONTINUE
RETURN
C
20 FORMAT (' STEP SIZE LIMITED AT Z=',F10.5)
END
C
C
C
SUBROUTINE GF (ZK,CO,SI)
C
C GF COMPUTES THE INTEGRAND EXP(JKR)/(KR) FOR NUMERICAL INTEGRATION.
C
REAL*8 ZK,CO,SI,RK,RKS,RKB2
COMMON /TMI/ ZPK,RKB2,IJ
C**
ZDK=ZK-ZPK
RK=DSQRT(RKB2+ZDK*ZDK)
SI=DSIN(RK)/RK
IF (IJ) 1,2,1
1 CO=DCOS(RK)/RK
RETURN
2 IF (RK.LT..2D0) GO TO 3
CO=(DCOS(RK)-1.)/RK
RETURN
3 RKS=RK*RK
CO=((-1.38888889D-3*RKS+4.16666667D-2)*RKS-.5)*RK
RETURN
END
C
C
C
SUBROUTINE TEST (F1R,F2R,TR,F1I,F2I,TI,DMIN)
C
C TEST FOR CONVERGENCE IN NUMERICAL INTEGRATION
C
REAL*8 F1R,F2R,TR,F1I,F2I,TI,DMIN,DEN
DEN=ABS(F2R)
TR=ABS(F2I)
IF (DEN.LT.TR) DEN=TR
IF (DEN.LT.DMIN) DEN=DMIN
IF (DEN.LT.1.D-37) GO TO 1
TR=ABS((F1R-F2R)/DEN)
TI=ABS((F1I-F2I)/DEN)
RETURN
1 TR=0.
TI=0.
RETURN
END
C
C
C
SUBROUTINE SFLDS (T,E)
C
C SFLDX RETURNS THE FIELD DUE TO GROUND FOR A CURRENT ELEMENT ON
C THE SOURCE SEGMENT AT T RELATIVE TO THE SEGMENT CENTER.
C
REAL*8 T,PI,TP,THET,POT,RK,R1,R2,ZMH,ZPH,SFAC
COMPLEX*16 E,U,U2,XX1,XX2,ERV,EZV,ERH,EZH,EPH
COMPLEX*16 T1,EXK,EYK,EZK,EXS,EYS,EZS,EXC,EYC,EZC,
1 ZRATI,ZRATI2,FRATI,ER,ET,HRV,HZV,HRH
INTEGER*4 IND1,IND2
COMMON/DATAJ/S,B,XJ,YJ,ZJ,CABJ,SABJ,SALPJ,EXK,EYK,EZK,EXS,EYS,
1 EZS,EXC,EYC,EZC,RKH,IND1,IND2,IPGND,IEXK
COMMON /INCOM/ XO,YO,ZO,SN,XSN,YSN,ISNOR
COMMON /GWAV/ U,U2,XX1,XX2,R1,R2,ZMH,ZPH
COMMON /GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
1 IFAR,IPERF,T1,T2
DIMENSION E(9)
DATA PI/3.141592654D0/,TP/6.283185308D0/,POT/1.570796327D0/
C**
C E WRITE(*,*) ' SFLDS: START'
C**
XT=XJ+T*CABJ
YT=YJ+T*SABJ
ZT=ZJ+T*SALPJ
RHX=XO-XT
RHY=YO-YT
RHS=RHX*RHX+RHY*RHY
C RHO=DSQRT(RHS)
RHO=SQRT(RHS)
IF (RHO.GT.0.) GO TO 1
RHX=1.
RHY=0.
PHX=0.
PHY=1.
GO TO 2
1 RHX=RHX/RHO
RHY=RHY/RHO
PHX=-RHY
PHY=RHX
2 CPH=RHX*XSN+RHY*YSN
SPH=RHY*XSN-RHX*YSN
IF (ABS(CPH).LT.1.E-10) CPH=0.
IF (ABS(SPH).LT.1.E-10) SPH=0.
ZPH=ZO+ZT
ZPHS=ZPH*ZPH
R2S=RHS+ZPHS
C R2=DSQRT(R2S)
R2=SQRT(R2S)
RK=R2*TP
XX2=CMPLX(DCOS(RK),-DSIN(RK))
IF (ISNOR.EQ.1) GO TO 3
C
C USE NORTON APPROXIMATION FOR FIELD DUE TO GROUND. CURRENT IS
C LUMPED AT SEGMENT CENTER WITH CURRENT MOMENT FOR CONSTANT, SINE,
C OR COSINE DISTRIBUTION.
C
ZMH=1.
R1=1.
XX1=0.
C**
CALL GWAVE (ERV,EZV,ERH,EZH,EPH)
C**
ET=-(0.,4.77134)*FRATI*XX2/(R2S*R2)
ER=2.*ET*CMPLX(1.,RK)
ET=ET*CMPLX(1.-RK*RK,RK)
HRV=(ER+ET)*RHO*ZPH/R2S
HZV=(ZPHS*ER-RHS*ET)/R2S
HRH=(RHS*ER-ZPHS*ET)/R2S
ERV=ERV-HRV
EZV=EZV-HZV
ERH=ERH+HRH
EZH=EZH+HRV
EPH=EPH+ET
ERV=ERV*SALPJ
EZV=EZV*SALPJ
ERH=ERH*SN*CPH
EZH=EZH*SN*CPH
EPH=EPH*SN*SPH
ERH=ERV+ERH
E(1)=(ERH*RHX+EPH*PHX)*S
E(2)=(ERH*RHY+EPH*PHY)*S
E(3)=(EZV+EZH)*S
E(4)=0.
E(5)=0.
E(6)=0.
SFAC=PI*S
SFAC=DSIN(SFAC)/SFAC
E(7)=E(1)*SFAC
E(8)=E(2)*SFAC
E(9)=E(3)*SFAC
C**
C E WRITE(*,*) ' SFLDS: RETURN LINE 91'
C**
RETURN
C
C INTERPOLATE IN SOMMERFELD FIELD TABLES
C
3 IF (RHO.LT.1.E-12) GO TO 4
THET=DATAN(ZPH/RHO)
GO TO 5
4 THET=POT
5 CALL INTRP (R2,THET,ERV,EZV,ERH,EPH)
C COMBINE VERTICAL AND HORIZONTAL COMPONENTS AND CONVERT TO X,Y,Z
C COMPONENTS. MULTIPLY BY EXP(-JKR)/R.
XX2=XX2/R2
SFAC=SN*CPH
ERH=XX2*(SALPJ*ERV+SFAC*ERH)
EZH=XX2*(SALPJ*EZV-SFAC*ERV)
EPH=SN*SPH*XX2*EPH
C X,Y,Z FIELDS FOR CONSTANT CURRENT
E(1)=ERH*RHX+EPH*PHX
E(2)=ERH*RHY+EPH*PHY
E(3)=EZH
RK=TP*T
C X,Y,Z FIELDS FOR SINE CURRENT
SFAC=DSIN(RK)
E(4)=E(1)*SFAC
E(5)=E(2)*SFAC
E(6)=E(3)*SFAC
C X,Y,Z FIELDS FOR COSINE CURRENT
SFAC=DCOS(RK)
E(7)=E(1)*SFAC
E(8)=E(2)*SFAC
E(9)=E(3)*SFAC
C**
C E WRITE(*,*) ' SFLDS: RETURN LINE 125'
C**
RETURN
END